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ABSTRACT 

The irradiation of protoplanetary discs by central stars is the main heating mech- 
anism for discs, resulting in their flared geometric structure. In a series of papers, 
we investigate the deep links between 2D self-consistent disc structure and planetary 
migration in irradiated discs, focusing particularly on those around M stars. In this 
first paper, we analyse the thermal structure of discs that are irradiated by an M 
star by solving the radiative transfer equation by means of a Monte Carlo code. Our 
simulations of irradiated hydrostatic discs are realistic and self-consistent in that they 
include dust settling with multiple grain sizes (N=15), the gravitational force of an 
embedded planet on the disc, and the presence of a dead zone (a region with very low 
levels of turbulence) within it. We show that dust settling drives the temperature of 
the mid-plane from an r~^/^ distribution (well mixed dust models) toward an r~^/^. 
The dead zone, meanwhile, leaves a dusty wall at its outer edge because dust settling 
in this region is enhanced compared to the active turbulent disc at larger disc radii. 
The disc heating produced by this irradiated wall provides a positive gradient region 
of the temperature in the dead zone in front of the wall. This is crucially important 
for slowing planetary migration because Lindblad torques are inversely proportional 
to the disc temperature. Furthermore, we show that low turbulence of the dead zone is 
self-consistently induced by dust settling, resulting in the Kelvin- Hclmholtz instability 
(KHI) . We show that the strength of turbulence arising from the KHI in the dead zone 
is a = 10"^ 

Key words: accretion, accretion discs - radiative transfer- turbulence - meth- 
odsmumerical - planetary systems:protoplanetary discs 



1 INTRODUCTION 

The discovery of a s as giant planet with a s mall orbital ra- 
dius around 51 Peg l|Mavor fc Quelo3ll995l ) motivated the 
proposal that such planets m ust form far from their cen- 
tral star, and then migrate l|Goldreich fc Tremaind 1 19791 . 

inward to their observed orbits. The total number 
of detected exoplanets has increased to more than 3000 
providing abundant supportive evidence for this idea. In 
addition to Hot Jupiters, massive Earth-like planets (1- 
10 M0) called Super-Earths, have re cently been discovered 
around M stars (|Beaulieu et al.ir2006l ). The detected number 
of such low mass exoplanets has also been increasing. These 
very different planetary populations both require migration 
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as an explanation of their observed mass - period relation 
dUdrv fc Santosll2007l '). 

It is well established that the thermal structure of proto- 
planetary discs plays a central role in determining planetary 
migration. This is because Lindblad resonances which con- 
trol the strength of the disc-planet inter action are strongly 
influenced b y the thermal gas pressure (|Artvmowic j Il993l : 
IWardlll997l ). The dominant heat source for discs is stel- 
lar irradiation, although viscous heating dominates within 
the^ central 1 au for classical T Tauri star ( CTTS) systems 
(i Chiang fc Goldreich 1997l, hereafter CG97; iD'Alessio 
ll'g'gsl ). The irradiation of the disc under the assumption of 
vertical hydrostatic equilibri um results in a flared disc shape 
teenvon fc HartmannI 0*9871 ) that reproduces the observed 
infrared excess in the spectral energy distributions (SEDs). 

Dust in protopl anetary discs is the m ain absorber of 
stellar radiation fe.g.. lDullemond et al.|[2009l ). An important 
ingredient in how discs are heated is the grain size distribu- 
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tion. iMathis et all (|l977l ) found that the size distribution of 
dust in the interstellar medium (ISM) is well represented by 
a power-law, the so-called MRN distribution. A related, but 
somewhat shallowe r power-law is likely to be applicable to 
dust in discs (e.g., iD'Alessio et al.ir200ll . hereafter DCHOl) 
due to grain growth, as discussed below. The composition 
of dust in discs is different from that in th e ISM mainly 
with respect to the cont ribution of carbon (Draine. fc Led 
11984 iPollack et al.|[l99l . hereafter P94). For discs, carbon 
is included as organics while they are graphite in the ISM. 
Recently, polycyclic aromatic hydrocarbons (PAHs) have re- 
ceived a lot of attention, and found that they pl ay an im- 
portant role in the SEDs (e.g.. [Praine fc Lill2007h although 
their formation mechanisms remain to be investigated. 

Another important aspect of dust in discs is its time 
evolution, that is, grain growth and dust settling. The 
sizes of dust grains gro w with time, resultin g in ever 
increasi ng dust settling (iDubrulle et al.l 1 19951. hereafter 
DMS 95; ISchrapler fc Henning|l2004l : iFromang fc Papaloizoul 
l2006h . More millimeter and less far-infrared emission 
in the observed SEDs are a good indication of grain 
growth and the resultant dust settling. The data have 
been re cently confronted b y comparisons with theoretical 
models (IChiang et al ] I2OOII : iDuUemond fc DominikI l2004al : 
ID'Alessio et al. I2OO6I ). Grain growth and dust settling are 



robust features of discs aroun d various ma sses of stars 
from Herbig Ae/Be (lAcke et al.|[2004.') . CTTSs (iFurlan et al 
200E ) to brown dwarfs (e.g., Apai et aLll2005l : Scholz et al 



20071 . hereafter S07). Grain growth is a lso thought to be the 



trigger of planetary format ion in discs (|Youdin fc GoodmanI 
120051 : 1 Johansen et al.llioOTl 'l. 

The structure of protoplanetary discs is far more inter- 
esting than just pure power-law behaviour, however. A ro- 
bust feature in discs is the presence of regions of very low am- 
plitude turbulent regions, called dead zones (jCammic 1996). 
Turbulence is most likely the outcome of the magn etoro- 
tational instability (MRI) (jBalbus fc HawlevllT991al lbh. The 
MRI requires good coupling between magnetic fields and the 
weakly conducting ionized molecules produced by the X-rays 
from the central star and cosmic-rays. The high density mid- 
plane region within 0.01 - 10 au of a star has difficulty in 
being ionized, resulting in a dead zone. Extensive studies 
so far have shown that dead zones are i nitially extended 
over roughly 10 au from the central star ( jSano ct al. 2000': 
llnutsuka fc Sandl2005l : iMatsumura fc Pudr itz 2003, 2006). 
They subsequently shrin k in time due to viscous evolution 
l|Matsumura et aUlioogh . 

Dead zones play a significant role in dust settling and 
the thermal structure of discs. This is because dust settling 
is determined by the balance between the gravity and tur- 
bulence (e.g., DMS95). Dust settling is enhanced in dead 
zones because turbulence is significantly reduced there. As 
we show below, this leaves a dusty wall at the boundary be- 
tween active and dead regions for turbulence. We will show 
that the direct irradiation of this high scale height of dust 
in the active region in turn creates a positive gradient of 
the disc temperature. We also demonstrate that dust set- 
tling in the dead zone drives a low level of turbulence by the 
Kelvin-Helmholtz instability with a turbulent amplitude of 
Q ~ 10-^ 

The presence of planets also distorts the di s tribu t ion of 
gas and dust in discs. Ijang-Condell fc Sasselovl (|2003l . 12004 



hereafter JS03, JS04, respectively) investigated the effect of 
the gravitational force of a planet on the t herma l structure 
of discs by using models of I Gal vet et al.l (|l99l|) . They as- 
sumed that stellar irradiation and viscous heating were the 
main heat sources where dust was assumed to be well-mixed 
with the gas. They found that the gravity of the planet pro- 
duces a self-shadowing (on dayside of planet) and a result- 
ing illumination region (on nightside of planet) around the 
planet. This causes temperature variations compared with 
the unperturbed case. The maximum variation is about 30 
per cent for planets which have the threshold masses to open 
up a gap. Recently, Ijang-Condelll (|2008l . hereafter JOS) im- 
proved the calculations of JS03, JS04 by considering the 
full 3D discs and including self-consistent treatments for the 
density and temperature structures, and found that the self- 
consistency is very important to determine the temperature 
structure accurately. This is because the temperature is very 
sensitive to the density distribution of dust. 

Thus, the dust distribution is crucial for the ther- 
mal structure of discs. A major consequence is its 
role in controlling planetary migration. Analytical and 
numerical simulations of this process so far assume 
discs to have is othermal or a sim p le power-law tem- 
perature profile (iTanaka et all |2002| : iNelson et all I2OO0I : 
iD'Angelo eral]|2003l ). More recently, non-isothermal discs 
have been treated in numerical sim ulations although they 
cannot include st e llar radiation (ID'Angelo et al.l l2003l : 



iKlahr fc Klevll2006l: [Paardekooper fc Mellemal I2OO6I ) . Ana- 
Ivticallv. lMenou fc GoodmanI ||2004[) have taken into account 
the radiation of stars, and Jang-Condell fc Sasselovl (|2005| ) 
have also included viscous heating and the gravity of plan- 
ets as well as stellar radiation although their disc models 
assume dust to be well-mixed with the gas, which is still 
far from realistic disc models. We have recently shown that 
dust settling has a very import ant effect on planetary mi- 
gration in a companion paper (iHasegawa fc Pudritz! f2009l . 
submitted) . 

In this paper, we present detailed radiative transfer cal- 
culations of the effects of dust settling and dead zones on 
the thermal structure of discs. We performed numerical sim- 
ulations, solving the rad iative transfer equation by mea ns of 
a Monte Carlo method (jDullemond fc Dominil3l2004bl '). We 
included dust settling with multiple grains sizes, a planet, 
and a dead zone in our 2D disc models. We especially fo- 
cus on discs around M stars whose properties are similar 
to discs observed around brown dwarfs (S07). This is phys- 
ically interesting because the recent primary target for ob- 
servations is low mass planets like Super-Earths and the de- 
tection probability is the highest around low mass stars and 
also because viscous heating can be safely ignored. Also, the 
computational time for the Monte Carlo method increases 
for massive discs. Low mass discs around M dwarfs are there- 
fore an ideal target. We show that the presence of dead zones 
can have a very significant impact on the thermal structure 
of protoplanetary discs. 

Our plan of this paper is the following. In § [2l we de- 
scribe our numerical methods, disc models, and dust proper- 
ties. In § [3l we present our results and describe how each of 
component (dust settling, a planet, and a dead zone) and the 
combination of them affect the density and thermal struc- 
ture of discs. In § |4l we verify our findings by performing a 
parameter study for disc properties and discuss other effects 
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on the temperature, such as viscous heating and planetary 
accretion heating that are ignored in our calculations. In § 
O we discuss the Kelvin-Helmholtz instability that can be 
active in the dead zone. Finally, we present our conclusions 
in §11 



2 NUMERICAL METHODS AND DISC 
MODELS 

2.1 Numerical methods 

We performed numerical simulation s by using a 2D radia- 
tive tr ansfer code, called RADMC0 (|Dullemond fc DominikI 
l2004bl ). RADMC is a versatile and highly reliable code based 
on the Monte Carlo techn ique in which the fu ll radiative 
transfer equation is solved (jPascucci et al.ll2004l 'l. The basic 
principle of the Monte Carlo metho d is very simple, and is 
well discussed in the literature (e.g.. lWhitnev fc Wolfi|[2002l . 
references herein) . The Monte Carlo method traces photons 
as they scatter and perform random walks through discs. 
The implementation of absorption processes is crucial for 
determining the the rmal structure of dis c s. In this code, an 
improved version of lBiork man fc Wood! (|200ll ) is adopted. 
Their method assumes discs to be in the local thermody- 
namic equilibrium, so that photons absorbed by the dust 
grains are forced to re-emit immediately. Furthermore, the 
re-emission of photons is adjusted so that they have the cor- 
rected temperature distribution. This makes it possible to 
avoid time-consuming iterative calculations. 

RADMC uses spherical coordinate systems in order to 
incre ase the efficiency of followi ng random walks of pho- 
tons (|Dullemond fc Turolla|[2000l ). However, cylindrical co- 
ordinate systems are generall y used in a nalytical theories 
of planetary migration (e.g., WardI [iQOTI ) . In order to es- 
timate the migration time accurately, we performed a co- 
ordinate transformation from spherical to cylindrical ones 
(see Appendix A). This is non-trivial because the path of 
photons projected into 2D planes are curved due to this 
projection effect, and because it depend s on projected co- 
ordinate systems l|van Noort et al.]|2002l ). Our implementa- 
tion of this coord inate transform was tested against a stan- 
dard benchmark (I Pascucci et al.ir2004 ). and proves success- 
ful (,Hasegawa.200i v 



2.2 Stellar & disc models 

We adopt a disc model that can reproduce the observed 
SEDs for M stars (S07). These models assume that the den- 
sity distribution of discs (gas + dust) can be represented 
by 



P = 



2Tvh 



exp 



2h^ J ' 



(2.1) 



of grain size distributions: one of which has large ho for a 
ISM-like grain size distribution (identical to that of gas), 
the other of which has small ho for a size distribution ex- 
tending to larger grain sizes, we use only large ho since we 
use much more detailed dust settling models, as discussed 
below. This large scale height is employed for gas and all 
sizes of dust when dust is assumed to be well-mixed with 
the gas, while it is used for gas and only the smallest dust 
grain when dust settling is included in disc models. For the 
star, we assume that it is described by blackbody radiation 
with the characteristic temperature T, . 

Table [T] summarises the stellar and disc parameters we 
use (S07). Note that we increase the total disc mass by a 
factor of 10 relative to that of S07 follow ing the literature 
l|Kokubo fc Idalll998l : fAlibert et ai]|2005h . This increment 
allows discs to form gas giants which have been mainly ob- 
served around ma ny intermediate massive stars like our Sun 
l|Udrvet al.ll2009l ). We run simulations only of the inner re- 
gion of the discs, which is typically half the standard size, 
in order to shorten computational time. This requires the 
disc mass to be adjusted, so that the surface density is kept 
10 times more massive compared with that of S07. We call 
this configuration our fiducial disc model. In § |4l we per- 
form parameter studies showing that our findings are model- 
independent. 



2.3 Heat sources 

We assume the main heat source for discs to be stellar irra- 
diation and neglect any other possible heat sources such as 
viscous heating. 

Viscous heating arises from the accretion processes of 
discs by stars. The temperature due to viscous heating under 
the assumption of gray atmosphere is 



3F^ 

4(73 



2 



(2.2) 



where as is the Stefan-Boltzmann constant, is the optical 
depth of the thermal emission of dust, the viscous flux Fv is 



3GM,Ma 
47rr3 



1/2 



(2.3) 



with E oc r^^ and h — ho{r/Rt,)^ (see Table[l]for the mean- 
ing of symbols). By adjusting ho and /3 with the flxed total 
disc mass and size, they reproduce the SEDs of very low 
mass stars such as M and brown dwarfs observed by Spitzer. 
Although they also took into account the effect of dust set- 
tling by considering two kind scale heights h for two kinds 

Our preliminary simulations showed that viscous heating is 
See the website: [http: / /www.mpia-hd.mpg.de/homes/duIIemon/r admc^miiBafatiatLly within 0.1 au for discs around M stars. Thus, 



and Mg is the accretion rate of discs by the star l|Pringlel 
Il98ll . also see Table [T}. For discs around the CTTSs, viscous 
heating dominates the heating by the star within about 1 
au (e.g., JSQ4). This is because the temperature due to stel- 
lar irradiation is roughly r~^/^ since the inverse square law 
while the temperature due to viscous heating is about r""^''* 
(see equation (|2.2p and (|2.3p ). Thus, the viscous tempera- 
ture decreases rapidly. 

For discs around M stars, however, the turnover point 
is at a much smaller disc radius than the case of CTTSs 
because the accretion rate Ma of M stars is a few orders 
of magnitude smaller than that of CTTSs. We a dopt Ma = 
10"" Mq yr-i for M stars (|Mohantv et al.ll2005h . Under the 
Eddington approximation, = 2/3 or higher corresponds 
to the optically thick regions. Thus, we consider = 1 
as the optically thick region although this choice does not 
change r„ very much since Tu oc t^'^ (see equation (|2.2p ). 
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Table 1. Summary of parameters & symbols 



Symbol 


Meaning 


Value 


M« 


Stellar mass 


0.1 Mq 


R* 


Stellar radius 


0.4Ro 


T, 


Stellar effective temperature 


2850 K 


Rout 


Outer disc radius 


50 au 


Rin 


Inner disc radius 


6R, 


To 


Characteristic disc radius 


R, 


ho 


Characteristic disc height 


0.02 To 




The exponent of disc height 


1.1 


S 


The surface density of (gas + dust) 


oc (r/ro)~l 


Ml 


The total disc mass (gas + dust) 


4.5 X 10~3 M0 




gas-to-dust ratio 


100 


a 


Parameter for turbulence (active/dead) 


10-2/10"^ 



M0 is solar mass, and Rq is solar radius. ^ is increased by a factor of 10 relative to S07. 



our neglect of viscous heating for low mass systems is valid, 
especially for M star systems. 



erally written as 



i(a) = no a^S{ai 



(2.4) 



2.4 Dust properties 

The thermal structure of discs is controlled by two main 
properties of dust, its composition and size distribution. For 
the composition of dust, we adopt the model of P94 in which 
real and imaginary refractive indices of dust are derived from 
available laboratory and astronomical data, theory, and the 
chemical composition of the primitive bodies in the solar 
system. Tabled summarises the composition, abundance, 
and density of dust. The opacity of dust is calculated based 
on the Mie theory using these refractive indices. 

For the composition, we assume that water ice exists 
over the entire extent of discs. This is because DCHOl found 
that mass absorption coefficients with and without the wa- 
ter ice are very similar (see their fig. 1), and because JS04 
found that the resultant density structures are also simi- 
lar (see their fig. 3). Thus, our assumption is appropriate. 
Note that the so-called ice line that defines the region in 
which wat e r ice can exist is important for planet formation. 
Ilda fc LinI l|2008l ) found that planet formation is accelerated 
at disc radii beyond the ice line, due to the formation of wa- 
ter ice and its retent ion arising from the pressure maximum 
iKretke fc Linll2007l ). This retention is likely to arise at the 
inner edge of the dead zones because the pressure maxima 
require a positive surface density gradient while we focus on 
the outer edge of the dead zones in this paper. 

For the size distribution of dust, we adopt discretised 
version of a power -law similar to the MRN distribution 
jMathis et al.l [l977l ') because dust settling causes different 
sizes of du st to h ave d ifferent scale heights. Furthermore, as 
shown by IWolj (|2003i ). the treatment of a real size distri- 
bution is important for calculating the dust temperature 
accurately compared with the approximate treatment in 
which the mean opacity for the ensemble size distribution 
is adopted. In contrast, we adopt the mean opacity for the 
composition. Thus the size distribution function n(a) is gen- 
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where no is a constant independent of grain size a, 5{a) is 
the Dirac's delta function, is the grain size we pick up, 
amin is the minimum size of dust, and amax is the maximum 
size of dust. We set amin = 0.01 (i.m and amax = 1000 (i.m 
(|Chiang et al.ll2001^ . The total dust mass mdust over the 
size distribution above becomes 



da— a pn(a) oc a 



3+p 



(2.5) 



where dust is assumed to be spherical and p is the mean den- 
sity for the ensemble of dust in the composition. In order for 
mdust to be an increasing function of a, we set p = —2.5. 
This shallower slope does not affect mass absorption coeffi- 
cients very much and is likely to be preferred in discs in order 
to take into account grain growth (DCIfOl). The intermedi- 
ate sizes of dust are logarithmically assigned between amin 
and amax- Fig. [T] depicts the maximum difference defined as 

'T{xo,N)-T{xo,N = 15) 



The maximum Diflt [%] = Max 



T{xo,N = 15) 



X 100 



(2.6) 

where T is the disc temperature, xo is the fixed position, TV 
is the total number of sampled grain sizes. We define the 
disc temperature T as below 

fa'pn{a)T{a) 



T = 



(2.7) 



where T{a) is the dust temperature for size a. The tempera- 
ture at A'^ = 15 is used as the reference value. The top panel 
shows the radial differences at the mid-plane and z = 0.1 au 
{— Xo) by the solid and dotted lines, respectively, and the 
bottom panel shows that the vertical differences at r = 1 and 
5 au {— Xo) by the solid and dotted lines, respectively. Both 
panels show the sign of convergence as more sizes of dust 
are included and that the maximum differences decrease to 
about or less 10 per cent. Our numerical convergence stud- 
ies show that 15 sizes of dust are sufficient and we adopt 
this in all models shown the paper. We emphasise that our 
calculations using a real size distribution enables one to pro- 
vide more accurate dust temperatures than using the mean 
iesjZip^iMtijtifo/otlaieifeiBsdiafali; of many dust grains done by S07. 
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Table 2. The composition and related quantities of dust 



Composition 


Abundance 


Density [g cm ^] 


Olivine 


2.6 


3.49 


Orthopyroxene 


0.8 


3.4 


Iron 


0.1 


7.87 


Troilite 


0.8 


4.83 


Organics 


4.1 


1.5 


Water ice 


5.6 


0.92 





Figure 1. The maximum temperature difference as a function 
of the total number of sizes of dust grains. Top: the radial tem- 
perature difference. The solid line denotes the difference at the 
mid-plane. The dotted line denotes the difference at z = 0.1 au 
from the mid-plane. Bottom; the vertical temperature difference. 
The solid line denotes the difference at r = 1 au from the central 
star. The dotted line denotes the difference at r = 5 au. For both 
cases, the temperature with 15 sizes of dust included is used as 
the reference. Both panels show convergence as more sizes of dust 
are added. 



2.5 Dust settling 

The motion of dust is diflFerent from that of gas. This is 
because dust is in Keplerian motion while gas is in slightly 
sub-Keplerian since it is affected by its thermal pressure. 
This difference implie s collisions or friction between them 
l|Weidenschilling|ll977^ ■ resulting in an identical density dis- 
tribution for them when the friction is efficient. Dust in pro- 
toplanetary discs is generally in the so-called Epstein regime 
wherein Tf /QJ,], < 1, where the friction timescale r/ is 



El — 

Pa Cs 



ps is the material bulk density of dust, pg is the gas density, 
a is the radius of dust grains, and Cs is the sound speed. 
In this regime, the coUisional efficiency (oc 1/t/) is mainly 
dependent on a. Thus, larger dust grains feel the friction 
from the gas less than smaller ones, resulting in decoupling 



from the gas. For laminar flows, the dust subsequently sinks 
into the mid-plane. Consequently, different sizes of dust have 
different scale heights, all of which are less than that of gas. 
This settling motion is prevented in turbulent flows - such 
as the MRI-active outer regions of discs. Dust settling is 
diffusive in a turbulent environment. Thus, dust settling is 
controlled by the balance between the gravitational force of 
the star and the amplitude of the turbulence in the discs for 
a given size of dust. 

We adopt the analytical model of DMS95 for dust 
settling in w hich turbulence in discs is described by a- 
prescription (IShakura fc SunvaevI Il973l ). Also, its energy 
spectrum is assumed to be represented by the Kolmogorov 
type, that is, E{k) oc k , where k is the wavenumber, 
and 'yturb is dependent on the nature of the turbulence. Since 
Iturb generally has the value between 5/3 and 3 (DMS95), 
we set 7turi) = 2. This choice is not signiflcant for our results 
(see equation (|2.1Up ). The reduced scale height of dust hd 
due to the dust settling for its size a is 



hd(a) 



where 



H 



1 



1 + 7turb 



1/4 



2Tvpsa 



(2.9) 



(2.10) 



and h is the scale height of gas. Equation (|2.10|l was derived 
from the advection diffusion equation with the diffusion coef- 
ficient calculated from the properties of the turbulence, and 
confirmed by their numerical calculations. This approach 
has been recently verified by magneto hydrodynamical sim- 
ulations fPromang fc Papaloizoull2006l ). 

Thus, the density distribution of the size a of dust with 
settling is described by equation (|2.H) with the replacement 
of h with hd{a). 



2.6 The gravity of planets 

The gravitational force of a planet distorts the distribution 
of gas in the disc (JS03; JS04). This distortion is taken into 
account by extending the original vertical hydrostatic bal- 
ance equation to one that includes the gravitational field 
of the planet as well as the central star. Consequently, the 
density distribution of gas becomes 



p ^ po exp 



z 

'2h? 



(2.11) 

where /i = Mp/M*, Mp is the planetary mass, rp is the 
location of the planet from the host star, and the normal- 
ization constant po is chosen so that the density at 2 = 
corresponds to the unperturbed density (following JS03). 
Note that the hydrostatic assumption breaks down within 
the Hill radius ru ~ rp( Mp/M«)^^^, so that we interpola te 
the density in this region (|jang-Condell fc Sasselo^l2005h . 

Thus, for a disc model without dust settling, but with a 
planet, the density distribution of dust for any size a is iden- 
tical to that of gas (equation 1)2. 11|) ') while the distribution 
of the size a of dust for the case of dust settling is repre- 
sented by the same equation above with the replacement of 
h with hd{a) that is calculated by equation (|2.9p . Note that 
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the density distribution of gas with dust settling as well as 
a planet is still represented by equation (|2.1H) . 

2.7 Dead zones with dust settling 

Our disc models include dead zones. In this paper, we as- 
sume a characteristic radius for the dead zone to be 4 
au in our fiducial disc model. Many studies have shown 
that it is roughly 10 au in th e early phases of disc evo- 
lution. iMatsumura et all (|2009l ) have recently shown that 
dead zones evolve with time. The outer edge of the dead 
zone moves inward following the accretion of disc material 
onto the star. In order to take time-evolution into account, 
we perform parameter studies in § |31 confirming the same 
results. We adopt the Q-prescription for turbulence, so that 
whether discs are active or "dead", they have a characteris- 
tic value of a. For dead zones (r ^ 4 au), aoz ~ 10~^ while 
OLactive = 10~^ for active regions (r > 4 au). We show later 
in the paper that a = 10~^ in dead zones is self-consistent. 

Dead zones have higher degrees of dust settling than 
active regions. Since a in equation (|2.10|l takes two different 
values, depending on either active or dead zones, the density 
distribution of dust for this case is described by equation 
(|2.H) with the replacement of h with hd{a) and two values 
of a. If a planet is present, it is represented by equation 
(|2.11|) with the replacement of h with hd{a) and two values 
of a. The distribution of gas is denoted by equation (|2.1[) 
and (|2.11[) with no and a planet, respectively. 

2.8 Self-consistency & the disc temperature 

Self-consistency is important for the disc temperature as JOS 
discussed. In order to achieve it, iterative calculations are es- 
sential. Our Monte Carlo calculations are time-consuming. 
About a week is needed to run a single disc model. This 
is because the mid-plane region in a disc, which is impor- 
tant for torque calculations, is optically thick so that a large 
number of photons (~ 10*) is required to avoid the random 
noise. Therefore, we perform an iteration for each disc model 
above. The temperature of each size of dust is first calcu- 
lated for an initial distribution of dust. The distribution of 
gas is not important for the temperature calculations since 
its heat capacity is much less than that of dust. Assuming 
the temperature of gas to be equal to the mass-averaged 
dust temperature (see equation (|2.7p ). the new scale height 
of gas is calculated, allowing us to calculate the new scale 
height of dust. Consequently, the new density distributions 
of gas and dust are obtained. 



3 RESULTS 

3.1 Dust settling 

In this subsection, we focus on the effect of dust settling 
on the density and thermal structures of discs. Fig. [2] shows 
the results of our simulations of the thermal and density 
structures of each dust grain without and with of dust set- 
tling (the left and right columns, respectively). Grain size 
increases from 0.01 |i.m, 3.16 ^.m to 1 mm on the top to bot- 
tom panels. They show that the temperature of a dust grain 
decreases as its size grows and that dust settling provides 



a lower temperature in the mid-plane region. Also, larger 
grains settle into much fiatter layers because of their lower 
friction with the turbulent gas. 

Fig. [3] shows the results of our simulations of the disc's 
temperature and density structure of the dust component for 
disc models with and without of dust settling (the bottom 
and top panels, respectively). The figures present the mass- 
averaged disc temperature given by equation (|2.7[) . The com- 
mon feature of both cases is that the surface layer has a 
higher temperature than the disc mid-plane. This is because 
the surface layer is directly heated by the central star while 
the mid-plane is heate d mainly by the thermal em ission of 
dust as first noted by Keny on fc H artmannI (|l987l . also see 
CG97). This is confirmed by the structure of the tempera- 
ture contours. For the optically thin surface layer, the con- 
tours have spherical shapes, meaning that the stellar radia- 
tion dominates the thermal emission of dust. For the mid- 
plane region, the contours show straight lines that are par- 
allel to the density structure of the dust denoted by colors. 
In this regime, dust emission dominates stellar radiation. 
The physical picture that the mid-plane region is enclosed 
by the superheated layers is consistent with the analytical 
models proposed by CG97 although their disc models have 
only two characteristic temperatures, one of which is for the 
superheated layer, the other of which characterizes the mid- 
plane. The transition region between these two regimes is 
characterised by a steep vertical temperature gradient. 

Dust settling changes the temperature structure as well 
as the density structure of dust as is expected by equa- 
tion (|2.9[) (see Fig. |3] and \^ . The main effect of dust set- 
tling is that the surface region has a higher temperature 
while the mid-plane region has a lower temperature. This 
behaviour is explained by the fact that different sizes of 
dust have different scale heights. In the surface layer, larger 
dust grains are depleted, but smaller grains still exist. These 
are more important for absorbing stellar radiation. Thus, 
photons emitted by the star can be efficiently absorbed by 
the surface layer. On the other hand, the mid-plane re- 
gion has many sizes of dust due to settling, resulting in a 
higher density compared with the case of well- mixed case. 
Thus, the optical depth is increased and photons emitted 
from the star have more difficulty in passing through the re- 
gion. Eventually, the energy absorbed in the region becomes 
smaller, and the temperature is decreased. This difference 
of the temperature, in turn, affects the density distribution 
of the gas by reducing its scale height (the vertical lines 
in Fig. |4]). Thus, dust settling drives the structure of discs 
from fiared to fiatter shapes (see Fig. [3]). This is consis- 
tent with the results of time-de pendent dust settling models 
IIDullemond fc Dominik|[2"004al ') while we assumed dust set- 
tling to be in the steady state. 

Fig. [S] shows the temperature at the mid-plane as a 
function of disc radius. The top panel shows the well-mixed 
model in which the temperature profile is well represented 
by r~^^^. Thus, our resultant discs are less flared relative to 
CG97 because of the increase of disc masses. The bottom 
panel shows the case of dust settling in which the temper- 
ature profile is well represented by r~^^*. Interestingly, this 
power- law can be derived from flat disc models (CG97). This 
is another indicator of flatter disc shapes. 
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Figure 2. The density and temperature structures of dust with its size a without and with dust setthng on the left and right 
columns, respectively. The dust density is denoted by the colors. The scale is shown in the color bar in units of [g cm~^]. Top: 
a = 0.01 i^m. Middle: a = 3.16|jm. Bottom: a = 1mm. The dust temperature is a decreasing function of a. Dust settling provides a 
flatter shape for the distribution of larger grains and a lower temperature in the mid-plane region. 



3.2 The effects of planets 

In this subsection, we focus on the effect of the gravitational 
force of planets with various masses. The fiducial location 
of a planet is fixed as 6 an. Fig. [U] shows the results of our 
simulations of the thermal structure of discs and the den- 
sity structure of dust with planets. The masses of the plan- 
ets are 0.1, 1, 5, and 10 Af^ on the top to bottom panels, 
respectively. The main difference in the density structure 
of dust is the decrement of the density above the planet 
which is expected by equation (|2.11|l . This effect increases 
as the masses of planet increase. Another difference is the 
decrement of the temperature of the mid-plane region on the 
dayside of the planet as shown by the contour of 10 K. This 
effect is also enhanced with planetary masses since it arises 
from the gravitational force of the planet to compress the 
material above them, resulting in the increase of the den- 
sity of dust in the mid-plane region (also see Fig. IBlf) . This 
makes it more difficult for photons to penetrate and heat the 
region. The lower temperature was not found by JS03, JS04, 
and JOS because they included viscous heating although it 
is dominant only within a few au. 

Fig. [7] shows the temperature at the mid- plane for the 
case of 10 . As discussed above, the temperature on the 
dayside of the planet becomes lower than that of the unper- 
turbed case (also see Fig. [5]). Furthermore, the temperature 
on the nightside of the planet also decreases for the same 



reason. Interestingly, the temperature at the location of the 
planet has a sharp peak. This arises from the compression 
of material above the planet, so that mean free paths of 
photons become longer. Consequently, more photons with 
higher energy can reach to the mid-plane regions. Further 
detailed analyses and discussion are tossed into Appendix 
B. 

We now focus on the combined effect of dust settling 
and planets. Fig. |S] shows the results of our simulations of 
the thermal structure of discs and the density structure of 
dust with dust settling and planets. Planetary masses are 
changed from 0.1, 1, 5, to lOMgj on the top to bottom panels, 
respectively. The important consequence of the combination 
is a dip of the temperature contours at the location of the 
planets, which is clearly seen by the contour of 10 K for the 
case of the planet with 0.1 M^. The existence of this higher 
temperature region above the planets is explained by the 
same reason of the effect of dust settling. In other words, 
the lower density of dust produced by the combination of 
the planet and dust settling increases the mean free paths 
of photons. This allows more photons penetrate deeper into 
the mid-plane region. Consequently, the region above the 
planets has a higher temperature. 

Fig. in] shows the temperature distribution with 10 
planet at the mid-plane. A temperature peak at the co- 
orbital radii of the planet is enhanced due to the combined 
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Figure 3. The density structure of dust and the temperature 
structure of disc. The dust density is denoted by the colors. The 
scale is shown in the color bar in units of [g cm~'^]. The tem- 
perature is denoted by the contours in the unit of [Kelvin]. Top: 
dust is well-mixed with gas. Bottom: dust settling included. Dust 
settling makes the surface layer higher temperature and lower 
density and makes the mid-plane region lower temperature and 
higher density. It also makes the dust distribution geometrically 
flatter (also see Fig. |4]l . 



eflFect of the planet and dust settling, as discussed above. It 
is likely that this effect may drastically change the forma- 
tion of gas giants. This is because accretion of gas envelops 
surrounding planets is sensitive to the gas temperature. Fur- 
ther detailed analyses and discussion of the dust distribution 
around a planet are presented in Appendix B. 



3.3 Dead zones in models with and without 
planets 

The major effect on the thermal structure of a disc arises 
from a dead zone (which is 4 au in size in our model). Fig. 
1101 shows the results of our simulations of the thermal struc- 
ture of discs and the density structure of dust with dead 
zones and dust settling. The dead zone drastically changes 
the density distribution of dust because dust settling is en- 
hanced due to the low turbulence there (also see equation 
I2.9|l . Thus, larger dust grains accumulate at the mid-plane 
since the scale heights of such dust become much smaller. 
The accumulation of dust into the mid-plane is prevented 
in the active zone due to the stronger turbulence. Conse- 
quently, a distinct step in the scale height of dust - or "wall" 
- is left at the boundary between the active and dead zones. 
This dusty wall plays a significant role in the thermal struc- 
ture of discs since dust is the main absorber of stellar irra- 
diation. 

The temperature structure of discs is also substantially 
changed by the dead zone although the basic physics oc- 
curring within the dead zone is the same as that of dust 
settling. Thus, the surface layer has a higher temperature 
and the mid-plane region has a lower temperature. Also, 




8x10" 
6x10" 



E 

" 4x10" 



well — mixed 
dust settling 



well — mixed 
dust settling 



0.4 0.6 

z (AU) 



Figure 4. The temperature and dust density structure at r = 5 
au as a function of distance from the mid-plane. Top: the temper- 
ature behaviour. Bottom: the dust density profiles. The vertical 
lines denote the scale height of gas. In both panels, any solid lines 
denote the case of well-mixed dust, and any dotted lines denote 
the case of dust settling. The reduced scale height of gas for the 
dust settling case also indicates flatter disc shapes. 
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Figure 5. The temperature structure as a function of disc radius. 
Top: the well-mixed case. The solid line denotes the temperature 
of our calculations. The dotted line denotes the temperature of 
flared disc models (CG97). The dashed line denotes the best fit 
for our calculations. Bottom: the dust settling case. The solid line 
denotes the temperature of our calculations. The dotted line de- 
notes the temperature of flared disc models (CG97). The dashed 
line denotes the best flt for our calculations, which is also derived 
as flat disc models by CG97. 



these two regions are separated by a transition region that 
features a tight collection of the temperature contours. The 
transition region is more declined toward the mid-plane in 
the dead zone. This arises from the enhancement of dust 
settling within the dead zone. Thus, in the dead zone, the 
superheated layer expands and the mid-plane region shrinks. 
In contrast, the turbulence in the active region keeps dust 
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Figure 6. The density structure of dust and the temperature structure of disc with well-mixed dust (as Fig. |3]l. The blow-up 
versions near the planet are shown in the right column. The thick solid lines denote the Hill radius of the planets. Top; 0.1 A/^. 
Second: 1 Mgj. Third: 5 Mgj Bottom: 10 Mq. The presence of the planets produces a lower temperature of the mid-plane region in 
front of them due to the compression of material above the planets (see the contour of 10 K). 



aloft although dust settling is inevitable because of grain 
growth due to agglomeration. The dust wall formed at the 
outer edge of the dead zone is thus directly exposed to the 
stellar radiation. This is one of the most important findings 
in our study. Consequently, the strange shapes of tempera- 
ture contours sandwiched by the spherical shapes (the opti- 
cally thin, surface regions) and straight lines (the optically 
thick, mid-plane regions) are produced. 

We found another important effect of the dusty wall. 
Fig. [TT] shows the temperature of the mid-plane. Impor- 
tantly, the temperature in the vicinity of the dead zone in- 
creases toward the dusty wall. This unusual temperature be- 



haviour is also explained by the direct exposure of the wall to 
stellar radiation. The direct heating provides the wall with 
higher energy, resulting in a higher temperature (otherwise 
the region is heated only by the thermal emission of dust, 
which provides much lower energy). Since the region is opti- 
cally thick, the ener gy absorbed by the wall is tr ansferred by 
diffusion processes (|Hasegawa fc Pudrit j I2OO9I ') . Thus, the 
presence of the hot dusty wall provides a positive tempera- 
ture gradient. It is well-represented by an r^^^ temperature 
distribution which can be also derived from analytical calcu- 
lations. This temperature behaviour provides striking effects 
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Figure 8. The density structure of dust and the temperature structure of disc with dust settling as well as a planet (as Fig. [BJ. 
The combination of dust settling and planets produces higher temperature above the planets due to the longer mean free paths of 
photons in the region (see the contour of 10 K). 



on pla netary migration, which are addressed in a companion 
paper (jHaseeawa fc Pudrit dlioogl ). 

We now focus on the combined effects of the dead zone 
with dust settling and planets with 10 M^. We examine 
two extreme cases; one of which is that of a planet that is 
located outside the dead zone (rp = 6 au), and the other in 
which it is located ai rp — 2 au inside the dead zone. Fig. 
[12] shows our results for the temperature structure of discs 
and the density structure of dust. The common feature of the 
density structure is the compression of the density above the 
planet although the overall density structure is similar to the 
case only with the dead zone. The main difference between 
these cases is the effect of the planet on the mid-plane region 



of discs. For the planet outside the dead zone, the density 
is further distorted by the gravitational force of the planet 
while, for the inside case, the distortion by the planet is not 
clear. This is explained by the enhancement of dust settling 
in the dead zone. As discussed above, dust settles much more 
compactly into the mid-plane in the dead zone. The planet 
also compresses dust into the mid-plane. Thus, both the 
dead zone and the planet affect the density distribution in a 
similar way, although the former is global and the latter is 
local. Thus, the global effect of the dead zone dominates the 
local effect of the planet. The active zone, however, allows 
the perturbation of the planet to be shown clearly due to 
the stronger turbulence. Thus, the combination of the dead 
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Figure 7. The temperature structure with a 10 Mgj planet as 
a function of disc radius (as Fig. [Sjl. The location of the planet 
is indicated by the arrow. The planet produces a slightly lower 
temperature region in the vicinity of the planet and a sharp peak 
at its location due to the compression of material. 
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Figure 10. The density structure of dust and the temperature 
structure of disc with a dead zone and dust settling (as Fig. [Sjl. 
The size of the dead zone is 4 au. The rapid settling of the dust 
in the dead zone leaves a dusty wall at the inner edge of the 
active zone which is directly heated by the central star, resulting 
in the strange shapes of temperature contours sandwiched by the 
spherical and straight lines. 
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Figure 9. The temperature structure with a 10 A/q planet and 
dust settling as a function of distance from the central star (as 
Fig-EJ- The combination of the planet and dust settling enhances 
a sharp peak at the location of the planet. 
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Figure 11. The temperature structure with a dead zone and dust 
settling as a function of distance from the central star (as Fig. 
[Sjl. A positive temperature gradient is formed at the boundary 
between the active and dead zones due to the direct heating by 
the central star. This is explained by radiative diffusion of the 
energy absorbed by the dusty wall (the dashed-dot line). 



zone and the planet still allows dust to settle into the mid- 
plane further although the dead zone surpasses the effect of 
the planet. Furthermore, the formation of the dusty wall is 
never affected by the inclusion of a planet, provided that the 
planets is not placed in the immediate vicinity of the dead 
zone. 

The temperature structure is also changed by the com- 
bination of the two although it is roughly similar to the case 
with only the dead zone. The main difference arises when the 
planet is placed outside the dead zone. The compression of 
dust produces the self-shadowing region shown by the con- 
tours of 10 and 20 K of Fig. [12] (also see Fig. [lOj. There is 
no clear difference for the case of the planet inside the dead 
zone. This is explained by the fact that the dust distribu- 
tion with the planet is almost identical to that without the 
planet since the effect of the dead zone is too strong relative 
to the gravitational force of the planet. Thus, the physical 
picture that the dusty wall becomes thermally hot due to 
the direct heating of the star is true even if the perturbation 
of the planets is included. This is also confirmed by Fig. [13] 
which shows the temperature at the mid-plane. 



4 THE EFFECTS OF DISC MASS AND 
VISCOUS HEATING 

We have demonstrated a novel thermal effect due to the pres- 
ence of the dead zone, namely, the appearance of the dusty 
wall at the inner edge of the active zone and the resultant 
positive temperature gradient. In this section, we undertake 
parameter studies to confirm this finding for various disc pa- 
rameters. Also, we discuss the effects of other heat sources 
that have been neglected so far. 



4.1 Parameter study 

We have used disc models and parameters that reproduce 
the observed SEDs around M stars (S07) in which E oc 1/r. 
On the other hand, there is another famous disc model, so- 
called mini mum mass solar nebula (MMSN) model in which 
E oc r"^/^ (|Havashi|[l98ll '). MMSN models also succeed in 
reproducing the observed SEDs (e.g., CG97). Thus, proto- 
planetary discs are likely to have the range for the surface 
density from to r~^^'^. Thus, it is important to check our 
findings for MMSN models as well. Also, we have increased 
the disc mass by a factor of 10 which is favoured if massive 
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Figure 12. The density structure of dust and the temperature structure of disc with a dead zone, dust setthng, and the planet 
with 10 Mgj (as Fig. [H]!. Top: the orbital radii is 2 au. Bottom: the orbital radii is 6 au. Both panels show the compression of the 
density of dust above the planet. For the bottom panel, a self-shadowing is seen clearly by the contour of 10 K (compare with Fig. 
nop . The dusty wall is maintained even with the presence of the planet for both cases. 
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Figure 13. The temperature structure with a dead zone, dust 
settling, and the planet with 10 as a function of disc radius (as 
Fig. mil . The combined effects produce the positive temperature 
gradient as well as the sharp peak at the location of the planet 
at the position of the planet. 



plane ts are to be formed l|Kokubo fc "ldalll998l : lAiibert et al.l 
I2OO5I ). 

With this motivation, we performed a parameter study 
involving the slope of surface density and the disc mass. Ta- 



ble [3] summarises our runs (also see Table [T]). In the simula- 
tions, we include a 5 planet at 6 au to verify the effect of 
the gravitational force of the planet. Fig. ll4l shows the tem- 
perature of the mid-plane for the runs. Every case shows the 
positive temperature gradient although MMSN-4 provides a 
shallower slope. This shallower slope is explained by a higher 
temperature of the mid-plane for lower mass discs, resulting 
in positive temperature gradients that are somewhat dimin- 
ished by the original temperature of the mid-plane. Thus, 
the formation of the positive temperature gradient depends 
on the mass of discs, not the slope of the surface density. 
Furthermore, sharp peaks produced by the presence of the 
planet are also sustained. Thus, our findings of the dusty 
wall and the positive temperature gradient are robust al- 
though the slope of the temperature may be changed by the 
mass of discs. 

We now investigate the effect of the size of the dead 
zones. We have fixed the size of the dead zone as 4 au so far. 
Many previous studies, however, showed that the size de- 
pends on parameters they used although there is the agree- 
ment wit h each other that i ts size is roughly 10 au. Fur- 
thermore, iMatsumura et al] |200^ found that dead zones 
evolve with time and shrink following the accretion of gas 
onto the star. Thus, the size of the dead zones also has some 
uncertainties. In order to compare our findings for various 
sizes of the dead zones, we performed simulations with 1 
and 8 au size of the dead zones with dust settling. In the 
simulations, we use the parameters of Table [1] and do not 
include the planet since we have shown that the resultant 
effects of the dead zone and the planet are robust. Fig. 1151 
shows the temperature at the mid-plane. Both cases show 
the positive temperature gradient at the outer edge of the 
dead zones. Furthermore, the temperature behaviours are 
also well-represented by r^^^, meaning that the energy ab- 
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Table 3. The parameter study 



Run 


S 


Ma 




MMSN-3 


J.-3/2 


4.5x10- 


-3 


MMSN-4 


^-3/2 


4.5x10" 


-4 


S07-4 


r-1 


4.5x10" 


-4 



sorbed in the dusty wall due to the direct heating from the 
star is transferred by the diffusion processes. Thus, our find- 
ings are independent of the size of dead zones. 

We next turn to examining dead zones with a finite 
radial transition region between the dead and active zones. 
Although we have assumed so far that the active region sud- 
denly appears at various disc radii from the central star, the 
transition region will in general be diffused since the pres- 
ence of the dead zones activates various instabilities as dis- 
cussed below. Dead zones with a finite transition region can 
be generally written as 



a{r) = anz + 



O^active 



tanh 



Ar 



+ 1 



(4.1) 



where ro = 4 an is the fiducial size of the dead zone and 
Ar is the thickness of the transition region. Fig. \Wi shows 
our results - the thermal and density structures of dust and 
the temperature of the mid-plane as a function of disc ra- 
dius from the star (the left and right columns, respectively). 
The parameter Ar is changed from 1, 3, 5, and 10 h, where 
h is the scale height of gas, on the top to bottom panels, 
respectively. The positive temperature gradient appears in 
every case. Interestingly, the slope of the positive tempera- 
ture gradient is more or less represented by r^^^ although 
it becomes shallower as Ar increases. The difference of the 
slope is explained by the shape of the dusty wall. A gradual 
transition of the value of a for a finite value of Ar provides 
smoother density distribution of dust, resulting in a gentle 
temperature gradient. Furthermore, the location of the pos- 
itive temperature gradient is changed for these cases. The 
positive temperature gradient appears at the mid-plane re- 
gion closer to the star, as Ar increases. This is also under- 
stood by the structure of the dusty wall. For the dead zone 
with a thicker transition region, the dusty wall presents at 
smaller disc radii from the host star. 



4.2 The addition of viscous heating 

In our simulations, we have considered stellar irradiation 
as the main heat source for discs, so that other potential 
heat sources such as viscous heating and the accretion lu- 
minosity of planets have been neglected. This is because we 
have assumed discs to be passive following CG97. This as- 
sumption is likely to be appropriate since viscous heating is 
dominant only within 0.1 au for M star systems. Further- 
more, the resultant SEDs of pa ssive discs agree with the 
observed ones very weU (CG97; IChiang fc GoldreichI [19991 : 
IChiang et al. I I2OO1I ). On the other hand, the SEDs produced 
from the models including viscous heating as well as stel- 
lar irradiation are als o able to reproduce the observations 
l|D'Alessio et al.lll998l . Il999l ). Thus, it is worth analysing 
the effects of viscous heating on the dusty wall and the re- 
sultant temperature structures since the constraints arising 
from SEDs are necessary, but not sufficient conditions. 




Figure 14. The temperature structure with a dead zone, dust 
settUng, and a planet with 5 as a function of disc radius (as 
Fig.[Tlll. Top: MMSN-3. Middle: MMSN-4 Bottom: S07-4. Every 
panel shows a positive temperature gradient and a sharp peak at 
the location of the planet (rp = 6 au) although the slope of the 
temperature depends on the disc mass. 



The positive temperature gradient formed in front of 
the dust wall is maintained even if viscous heating is in- 
cluded. This is because the temperature around the dust 
wall determined by the direct exposure of stellar irradiation 
is much higher. Viscous heating cannot wash out this effect 
since the temperature due to viscous heating is dominant 
only within 0.1 au. 

Furthermore, the low turbulence within the dead zones 
implies that viscous heating is strongly reduced. Weaker vis- 
cous heating strongly supports our assumption of passive 
discs. Therefore, our assumption is valid throughout the disc 
beyond 0.1 au and our findings are robust features of the 
dead zone even if viscous heating is included. 

In addition to viscous heating, the presence of plan- 
ets provides another heat source due to accretion processes. 
During planet formation, protoplanets accrete gas, dust and 
planetesimals, resulting in the radiation of the extra gravi- 
tational energy. The temperature structure produced by the 
accretion luminosity of the planet is represented by r~^^^ 
due to the inverse square law. Thus, the temperature peak 
produced by the gravitational force of the planet becomes 
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Figure 16. The density structure of dust and the temperature structure of disc with the dead zone and dust setthng with a finite 
size of the transition region Ar (as Fig. |3)l and the temperature of the mid-lane as a function of disc radius (as Fig. on left and 
right columns, respectively. Top: Ar = Ih. Second: Ar = 3h. Third: Ar = 5h Bottom: Ar = lOh. Every case shows a dusty wall 
and the resultant positive temperature gradient although its slope becomes shallower and its position becomes closer to the central 
star, as Ar is increased. 



wider due to the accretion effect. Thus, although the accre- 
tion luminosity provides slight modification, the basic pic- 
ture of our results is still valid. Planetary migration, how- 
ever, may be affected by the accretion luminosity since it is 
determined by the quantities in the vicinity of planets. We 
leave this study in the future publication. 



5 INDUCING DEAD ZONE TURBULENCE BY 
DUST SETTLING 

The presence of a dead zone is known to activate various 
instabilities in protoplanetary discs. In this section, we focus 
on the Kelvin-Helmholtz instability (KHI) that may arise. 

The KHI is triggered by dust settling because the ver- 
tical sh«a£isjii^uced_b2_ths^ of dust den- 
sity (|Sekivalll998l : |Johansen et al.ll2006l '). The increase of the 
density of dust in the mid-plane causes the gas to be in Ke- 
plerian motion since the collisional interactions with dust 
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Figure 15. The temperature structure with a dead zone and dust 
settUng as a function of distance from the central star (as Fig. llll l. 
Top: the size of the dead zone is 1 au. Bottom: the size of the dead 
zone is 8 au. Both panels show a positive temperature gradient 
formed at the boundary between the active and dead zones. 



become more efficient while the gas above the mid-plane is 
still in sub-Keplerian due to lower density of dust. This ver- 
tical difference in the gas velocity induces a vertical shear 
flow, resulting in the KHI. Such regions can therefore be- 
come turbulent. 

This dust-induced turbulence is local and well- 
represented by the so-called Richardson number Ri, 



Ri 



gzdln{p + pd)/dz 
{duy/dzY 



(5.1) 



where is the vertical gravitational force, p is the gas den- 
sity, pd is the dust density and Uz is the vertical component 
of the gas velocity. The cla ssical theory finds tha t Ri < 1/4 
is the unstable condition (| Chandrasekharl Il96l[ ) although 
it is the necessary, but not sufficient condition. Recent nu- 
merical simulations h ave shown that the critical value is 1 
(|johansen et al.ll2006l ). Furthermore, the critical Ri allows 
one to estimate the critical scale height of dust hd,crit, 



hd 



' 2jad ' 



(5.2) 



where ^ = (h /r)(d\n p/d\ nr) and 7ad ~ 5/3 is the adiabatic 
index of gas (|Sekivalll99^ ). 

As discussed in § [2l dust settling is compensated by 
the turbulence of discs. In our disc models, the strength of 
turbulence is controlled by a which is a useful prescription 
for any global turbulence. Thus, we have treated the active 
and dead zones by adjusting a. In the active well coupled 
zone, the origin of the turbulence is obvious, that is, the 
MRI. The origin of the turbulence in the dead zone, however, 
is ambiguous since the MRI cannot be active due to the low 
ionization. Thus, we examine whether or not our assumption 
for the dead zones is consistent with the KHI. 



Fig.flTlshows the scale height of dust as a function of the 
distance from the star with a 4 au sized dead zone. Since the 
scale height of dust depends on the strength of turbulence a 
as well as its size a in our disc models (see equation (12. 9p ). it 
drastically changes at the boundary of the active and dead 
zones. Also, the critical scale heights of dust with Ri = 1/4 
and Ri — 1 are plotted. In the active zone, both of the 
critical scale heights of dust are much smaller than those 
of our scale heights of dust. This indicates that the MRI 
turbulence is so strong that the KHI cannot be active and 
is not therefore the primary driver of disc turbulence. 

The situation in the dead zone is significantly different. 
The critical scale height of dust with Ri — 1/4 occurs for 
dust grains with 1mm size. For the case of Ri = 1, the 
critical case is that for dust with 193 \Lm size. For these 
grains, the KHI becomes active. This implies that the origin 
of turbulence in the dead zones may be due to the KHI. 
Thus, it is wor th estimating the valu e of a arising from the 
KHI following [johansen et all (120061 ). From equation (|2.9|) . 
the ratio of the scale height of dust to that of gas becomes 



h 



(5.3) 



since H <ti 1. In addition, the scale height of dust estimated 
from the KHI is given by equation (|5.2p . Equating these two 
equations, the aKHi arising from the KHI is 



OiKHI 



(5.4) 



Fig. [TS] shows the value of okhi as a function of disc 
radius. As discussed above, the ukhi in the active zone 
is about two orders of magnitude smaller than that of the 
MRI. Within the dead zone, the Okhi lies in the range be- 
tween 10~^ and 10~*. This support s the value of a in the 
dead zone used in the literat ure (e.g.. llnutsuka fc Sanoll2005l : 
iMatsumura fc Pudritzl2006l '). The origin of turbulence in the 
dead zone is probably governed by the KHI, although the 
derived a is a kind of the maximum value because the recent 
study including the radial Keplerian shear has shown that 
the KHI can be active only when a mode of the instabilit y 
grows faster t han it is sheared out Jlshitsu fc" Sekivall2003l ). 
Furthermore, llshitsu fc Sekival (120031 7 have shown that al- 
though the equilibrium state of the dust distribution can 
be established, the critical conditions are different from the 
Richardson number because it does not include the effect of 
the Keplerian shear. Further study such as global simula- 
tions is needed. 



6 CONCLUSIONS 

In this paper, we have investigated the detailed thermal 
structure of protoplanetary discs around M stars by solv- 
ing the radiative transfer equation in 2D discs by means of 
the Monte Carlo method. In our simulations, stellar irradi- 
ation is absorbed by dust in discs and is their main heat 
source. Thus, the properties of dust such as its composition, 
size and density distributions are crucial in this study. For 
the composition, we have adopted the models that are de- 
rived from the solar abundance. For the size distribution, we 
have adopted the power-law distribution in which the size of 
dust is from 0.1 ^.m to 1 mm. Since we have used elaborated 
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Figure 17. The scale heights of dust as a function of distance 
from the central star with a 4 au sized dead zone. The scale heights 
of dust with various grain sizes in our disc models are shown. In 
addition, the critical scale heights with Ki = 1/4 and Ri = I 

are denoted by the solid and dotted lines, respectively. The KHI 
becomes active for dust grains with 1 mm and with 193 [xai when 
Ri = 1/4 and Ri = 1, respectively. 
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Figure 18. The values of a as a function of disc radius. The 
ctKHi estimated from the KHI with Ri = 1/4 and Ri = 1 are 
denoted by the solid and dotted lines, respectively. The value of 
a in our disc models is denoted by the dashed line. The value of 
a within the dead zone lies the range between 10"* and 10"^ for 
both values of Ri. This is consistent with the literature. 



models for the dust distribution, we had to sample the size 
of dust within this range. Our convergence study has shown 
that 15 sizes of dust chosen from this range are sufficient to 
describe the temperature of discs to an accuracy of better 
than 10 per cent. Our results of the dust temperatures are 
more accurate than those derived from the mean opacity for 
the ensemble of larger dust grains (e.g., S07). Furthermore, 
we have included dust settling, the gravitational force of the 
planet and the presence of dead zone. 

The most important findings of this study are the fol- 
lowings; 

(i) The temperature of the mid-plane is well-represented 
by r"^/* when dust settling is included. This power-law was 
derived by CG97 by assuming discs to be flat. Well mixed 

dust models result in an r^'^^^ temperature distribution. We 
have confirmed that dust settling makes discs geometrically 
flatter. 

(ii) The presence of the dead zone in discs implies the 



presence of a dust wall at the boundary of the active and 
dead zones. This wall will be directly exposed to stellar ir- 
radiation. Consequently, a positive temperature gradient is 
formed in the temperature profile of the mid-plane. This is 
explained by the higher energy at the wall that is propa- 
gated by radiative diffusion since the mid-plane is optically 
thick. We have found that even if the planets are present 
in front of the wall, this thermally hot dusty wall and the 
positive temperature gradient are maintained. 

(iii) The dead zones excite the Kelvin-Helmholtz instabil- 
ity because of enhanced dust settling within them. Conse- 
quently, we have demonstrated that a self-consistent value 
of Q ~ 10~^ in dead zones is expected. This result backs up 
the assumptions of many previous models for dead zones. 

Furthermore, we have confirmed earlier results in the 
literature that the surface layer has a higher temperature 
than the mid-plane region in any disc model. The surface 
layer is directly heated by the central star while the mid- 
plane region is only heated by the thermal emission of dust. 
Our simulations with dust settling have shown that this pro- 
cess leaves the surface layer hotter and the mid-plane region 
cooler. This arises from the deficit of larger dust grains in the 
surface layer and the increment of them in the mid-plane re- 
gion. This can be also understood by the mean free paths of 
photons that are determined by the density of dust. There- 
fore, dust settling makes the mean free paths in the surface 
layer longer and those in the mid-plane region shorter. Also, 
we have shown that the inclusion of the planets triggers the 
compression of the material above them, resulting in longer 
mean free paths of photons. The region above the planets 
has a higher temperature. 

We have also performed a parameter study for our disc 
models by varying the total disc masses and the surface 
density distributions. It is generally thought that the disc 
masses arc about ten times more massive than the observed 
ones in order for gas giants to be formed. Also the obser- 
vations have indicated that the surface density has some 
range from to r~^^^. We have confirmed that our find- 
ings are basically unaffected for such parameter spaces, al- 
though the lower disc masses result in a higher temperature 
at the mid-plane, causing the positive temperature gradient 
to be somewhat diminished. In addition, the formation of 
the positive temperature gradient is independent of the size 
of dead zones since the dusty wall is present for any size of 
the dead zone. Furthermore, the positive temperature gra- 
dient is sustained for dead zones with a finite transition for 
a. 

We have also examined other heat sources for discs such 
as viscous heating and the accretion luminosity of the plan- 
ets. We have estimated viscous heating is only important 
within about 0.1 au for M stars although it is known to be 
dominant within about 1 au for the classical T Tauri stars. 
Thus, we have investigated the detailed effects of dust set- 
tling, dead zones, and planets on the thermal structure of 
irradiated discs. 

In subsequent papers, we will discuss the effect of our 
results on planetary migration because Lindblad torques, 
the main driving force of type I migration, arc affected by 
the gas pressure. In addition, we will address the effects of 
the dead zones on the spectral energy distributions (SEDs) 
and images. It is well known that the presence of embedded 
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planets in discs provide some interesting observational sig- 
natures on th e SEDs and/or image s by forming a gap for 
massive ones (jVarniere et al.l l2006l ) a nd by di storting the 
density distribution for low mass ones (|jang-Condell. 20091 ). 
We anticipate that the presence of the dead zones will also 
have an interesting contribution to the SEDs arising from 
the dusty wall and resultant positive temperature gradient. 
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APPENDIX A: THE TRAJECTORY OF 
PHOTONS PROJECTED INTO 2D PLANES IN 
CYLINDRICAL COORDINATE SYSTEMS 

Here, we obtain solutions of the trajectory of photons pro- 
jected into 2D planes (r, z) in cylindrical coordinate sys- 
tems. First of all, we setup a global cylindrical (r, z, $) 
and local spherical [9, tj)) c oordinate system by following 
iDuUemond fc Turollal l|2000l ) . We assume discs to be axisym- 
metrical, so that any dependence on "I> is dropped out. The 
north pole is aligned with the z-axis of the global system. 
Consequently, 9 is measured from this axis toward the mid- 
plane and follows photons as they move upward or down- 
ward compared with the mid-plane, (j) is measured from the 
local X-axis which is parallel to the global x-axis. The (j) co- 
ordinate tracks photons as they move forward or backward 
compared with the local coordinate system. 

The radiative transfer equation in the cylindrical system 



the constants are 



dl^ r 5" . ,dlv , dl^ , x/1 ^ A*^ cos </) div 

— = Vl-P^sm(?!)— + — , 

as or oz r oq) 

(Al) 

where is the specific intensity at frequency v, ds is the dis- 
place ment of the path of a photon, and ^ — cos 9 (|Liu et al.l 
l2006t ). Thus, the variations of r, z, fi, and (j) along the path 
of the photon are 



dz 
ds 
djj. 
ds 



A*, 
0, 



ds r 
A set of solutions we found are 

= {l-^^)s'' + b\ 

z = /ioS + 20, 

MO = 



sm = 



(A2a) 
(A2b) 
(A2c) 

(A2d) 



(A3a) 
(A3b) 
(A3c) 

(A3d) 



where 6 is the impact parameter of the ray with respect to 
the origin, zo is the height from the mid-plane of the closest 
point to the north pole of the global system. 

We use the solutions to follow the path of photons. 
When a photon is located at P = (r^, zi) and its new direc- 
tion is determined as (/i, (ji) = (/li, i\>j) by a random number. 
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(A4a) 
(A4b) 



where fc, i, and j represent the grid number at which the 
photon is placed since although the Monte Carlo method is 
a non grid-based code, all physical quantities such as the 
density and temperature are defined on the grid. Thus, the 
next point it will reach is determined by the intersection 
between grid points and the solutions. 

Equations (|A3aP and (IA3b|) give us two and one values 
for s, respectively. Since possible grid points the photon can 
take are Vk-i, r^, and rt+i {rt-i < rt < r^+i), and 
Zl, and zi+i (zi-i < zi < zi+i ). Thus, the total number of 
possible values of s is 9 and are given by 



= ± 



fe2 



ZL 



ZO 



(A5a) 



(A5b) 



MO 

where n = 1, 2, 6, = {fc — 1, fc, fc -I- 1}, m = 7, 8, 9, and 
L = {I — 1,1, 1 + 1} . The actual path is determined so that 
the absolute value of s is a minimum and not the same as 
sp which is the path evaluated at P = {r^, zi). 

Thus, the motion of photons projected into 2D plane in 
cylindrical coordinate systems is determined by equations 
(|A5al) and (|A5b|) . 



APPENDIX B: TEMPERATURE AND DUST 
DENSITY DISTRIBUTION AROUND A 
PLANET 

Here, we carefully analyse and discuss the thermal and dust 
density structures of discs with a planet for the well mixed 
and dust settling cases. 



Bl The well-mixed case 

Fig. IB II shows the vertical temperature and density of dust, 
in a model with well mixed dust, in the vicinity of the 10 
planet (Solid lines). For comparison, the unperturbed case 
is denoted by the dotted lines. Quantities at — th are de- 
noted in red while those at rp-\-rH are in black. The presence 
of a planet increases the density at the mid-plane region due 
to the gravitational force of the planet. This results in the 
compression of dust (and gas) at rpiirn (compare solid lines 
with dotted lines). For the temperature structure, the inclu- 
sion of the planets causes the surface layer to be hotter at 
rp — Th than the case without the planet (see two red lines) . 
This is explained by similar effects of dust settling. The com- 
pression of material reduces the optical depth, resulting in a 
higher temperature. On the other hand, the surface layer at 
Tp + TH becomes slightly cooler by the planets (see two black 
lines). In other words, there is no self-illumination region in 
our simulations. For the mid-plane region, the temperature 
becomes cooler due to the gravitational force of the planets 
at Tp ± Th (compare solid lines with dotted lines) . This is a 
self-shadowing effect, as found by JS03, JS04, JOS. 

Thus, we could not find the self-illumination regions 
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Figure Bl. The temperature and dust density structure in a disc 
witli well mixed dust and a 10 planet as a function of distance 
from the mid-plane. Top: the temperature behaviours. Bottom: 
the dust density profiles. In both panels, solid lines denote the 
case with the planet, and dotted lines denote the case without the 
planet for comparison. The red lines denote quantities at Vp — rjj , 
and the black lines denote quantities at rp+rg. The perturbation 
due to the planet accumulates the material near the mid-plane 
region, and produces the self-shadowing region in the vicinity of 
the planet. 
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Figure B2. The temperature and density structure in a disc 
with dust settling and a 10 Afgj planet, as a function of distance 
from the mid-plane (as Fig. IBH . Higher dust density in the mid- 
plane region is due to the combined effect of the planet and dust 
settling. This results in the much smaller temperature difference 
between rp ± rj{ . 



clearly, differing from the findings of JS03, JS04 and JOS. 
We did find the self-shadowing regions, however, although 
its effects are relatively small. This is because of the incre- 
ment of disc masses and the resultant relatively less fiaring 
disc structure. JS04 found that the temperature variation 
due to the planets becomes smaller when the slope of the 
surface density is steeper because higher density prevents 
photons from penetrating it deeply. Also, JOS found that 
the self-shadowing and illumination regions are very sensi- 
tive to the incident angle of stellar radiation. The increase 
of mass allowing massive planets to be formed results in a 
relatively fiatter disc structure (see Fig. [5]). This prevents 
the regions perturbed by the planet from being exposed to 
stellar radiation. Thus, the results of our simulations are not 
precisely the same, although consistent with those of JS03, 
JS04 and JOS due to the increase of the disc mass. 



B2 The dust settling case 

Fig. IB2l shows the vertical temperature and density of dust, 
in a model with dust settling, in the vicinity of the planet 
with 10 M®. The main trend is the same as the case in- 
cluding the planets only. Higher density in the mid-plane 
regions results in much smaller temperature difference be- 
tween Vp ± th- For the surface region, the temperature at 
Vp — rn becomes higher while the one at rp + th is lower 
when the planets are included. 



This paper has been typeset from a TpjX/ file prepared 

by the author. 



